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Abstract We propose an efficient Markov Chain Monte Carlo method for sampling 
equilibrium distributions for stochastic lattice models, capable of handling correctly 
long and short-range particle interactions. The proposed method is a Metropolis- 
type algorithm with the proposal probability transition matrix based on the coarse- 
grained approximating measures introduced in [ 17 , 21 1. We prove that the proposed 
algorithm reduces the computational cost due to energy differences and has compa- 
rable mixing properties with the classical microscopic Metropolis algorithm, con- 
trolled by the level of coarsening and reconstruction procedure. The properties and 
effectiveness of the algorithm are demonstrated with an exactly solvable example 
of a one dimensional Ising-type model, comparing efficiency of the single spin-flip 
Metropolis dynamics and the proposed coupled Metropolis algorithm. 
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1 Introduction 



Microscopic, extended (many-particle) systems with complex interactions are ubiq- 
uitous in science and engineering applications in a variety of physical and chemi- 
cal systems, exhibiting rich mesoscopic morphologies. For example, nano-pattern 
formation via self-assembly, arises in surface processes e.g., in heteroepitaxy, in- 
duced by competing short and long-range interactions |6|. Other examples include 
macromolecular systems such as polymers, proteins and other soft matter systems, 
quantum dots and micromagnetic materials. Scientific computing for this class of 
systems can rely on molecular simulation methods such as Kinetic Monte Carlo 
(KMC) or Molecular Dynamics (MD), however their extensivity, their inherently 
complex interactions and stochastic nature, severely limit the spatio-temporal scales 
that can be addressed by these direct numerical simulation methods. 

One of our primary goals is to develop systematic mathematical and computa- 
tional strategies for the speed-up of microscopic simulation methods by developing 
coarse-grained (CG) approximations, thus reducing the extended system's degrees 
of freedom. To date coarse-graining methods have been a subject of intense focus, 
mainly outside mathematics and primarily in the physics, applied sciences and en- 
gineering literatures 1 10 17 25 28 30 1. The existing approaches can give unprece- 
dented speed-up to molecular simulations and can work well in certain parameter 
regimes, for instance, at high temperatures or low density. On the other hand, in 
many parameter regimes, important macroscopic properties may not be captured 
properly, e.g. p] [30][3T[ . Here we propose to, develop reliable CG algorithms for 
stochastic lattice systems with complex, and often competing particle interactions 
in equilibrium. Our proposed methodologies stem from the synergy of stochastic 
processes, statistical mechanics and statistics sampling methods. 

Monte Carlo algorithms provide a computational tool capable of estimating ob- 
servables defined on high-dimensional configuration spaces that are typical for mod- 
eling of complex interacting particle systems at or out of equihbrium. Markov Chain 
Monte Carlo (MCMC) simulation methods such as the Metropolis algorithm, were 
first proposed in 1953 by Metropolis and his coauthors |29J for the numerical cal- 
culation of the equation of state for a system of rigid spheres. It was generalized in 
1970 by Hastings fT?! and it is commonly referred to as the Metropolis-Hastings 
(MH) Monte Carlo method. This method belongs to the family of MCMC methods 
which generate ergodic Markovian chains with the stationary distribution being the 
desired sampled probability measure. Metropolis algorithm consists of two main in- 
gredients: (a) the probability transition kernel q\ the proposal, that generates trial 
states and (b) the acceptance probability a according to which the proposed trial 
is accepted or rejected. There are though some drawbacks of this method when ap- 
phed to large systems, such as a small acceptance probability a, that leads to costly 
calculations of a large number of samples that are discarded. A way to reduce these 
costs is to predict efficient proposal measures such that the computational cost of 
calculating a sample is lower and, if possible, increase the acceptance probability. 
Convergence and ergodicity properties of Metropolis type algorithms are studied 
extensively in a series of works ||7]|8]|32). The rate of convergence to stationarity 
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is strongly dependent on the proposal distribution and its relation to the stationary 
measure ( p2) ch. 7). A quantity that measures the speed of convergence in distri- 
bution to stationarity is the spectral gap. In order to improve an MCMC method one 
has to increase its spectral gap by smartly constructing a good proposal. 

In this work we propose the Coupled Coarse Graining Monte Carlo (Coupled 
CGMC) method, a new method of constructing efficient proposal measures based 
on coarse-graining properties of the sampling models. We prove that such approach 
is suitable for models that include both short and long-range interactions between 
particles. Long-range interactions are well-approximated by coarse graining tech- 
niques 1 17 18 20), and Coarse Graining Monte Carlo (CGMC) are adequate sim- 
ulation methods with an order of acceleration up to 0{q^) with q a parameter con- 
trolling the level of coarse graining | [T9l[2T| . Furthermore, models where only short- 
range interactions appear are inexpensive to simulate, for example with a single 
spin-flip Metropolis method. However, when both short and long-range interactions 
are present the classical MH algorithm becomes prohibitively expensive due to the 
high cost of calculating energy differences arising from the long-range interaction 
potential. In p3] we extend our framework for coupled CGMC to the dynamics 
case, developing Kinetic Monte Carlo algorithms based on coarse-level rates. 

Section|2]describes the classical Metropolis-Hastings algorithm and some known 
mathematical theory for convergence and the rate of convergence for MCMC meth- 
ods. In Section [3] we present the proposed Coupled CGMC method in a general 
framework describing its mathematical properties. We state the main theorem that 
compares the rate of convergence to equilibrium with the rate of the classical MH 
method. In Section |4] we describe stochastic lattice systems and the coarse-graining 
procedure in order to prepare for the application of the proposed method in Section|5] 
to a generic model of lattice systems in which both short and long-range interactions 
are present. 



2 MCMC methods 

Before describing the Metropolis-Hastings method we need to introduce some nec- 
essary definitions and theoretical facts. 

Let X„ be a Markov chain on space L with transition kernel J^. 

Definition 1 A transition kernel ,J/if has the stationary measure /i // 

JfTpL = 11. 

Definition 2 J(f is called reversible with respect to pi if 

[g, .Jeh)^,^ {.J^g, h)^ , for all g.heL^ill). 
where {g,h)^ = j^g{o)h{o)ii{da) and J(fg{o) = f^.^{a,da')g{a')ya € E. 
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A sufficient condition for jj. being a stationary measure of Jt^ is the, often easy 
to check, detailed balance(DB) condition. 

Definition 3 A Markov chain with transition kernel Ji(f satisfies the detailed bal- 
ance condition if there exists a function f satisfying 

J^{a,a')f{a)^Jf{a',a)f{a'). (1) 

Here we focus on the Metropolis-Hastings algorithm p2| . The algorithm gen- 
erates an ergodic Markov chain X„ in the state space E, with stationary measure 
H{da). Let /((t) be the probability density corresponding to the measure /i and 
Xq = Ob be arbitrary. The n-th iteration of the algorithm consists of the following 
steps 

Algorithm 1 (Metropolis-Hastings Algorithm) 



Given X„ = (7 

Step 1 Generate Y„ = a' ^ q{G'\G) 
Step 2 Accept-Reject 



where 



Y„ = o' with probability a((J,(y') 
= <7„ with probability 1 — CX{<J, o') 



, . /, ficy')qia',a)\ 
a{a, a') =mmi 1, - ) -f } 



We denote q{(y'\(y) the proposal probability transition kernel, and o;((T,(7') the 
probability of accepting the proposed state o' . The transition kernel associated to 
MH algorithm is 



■J€c{a, a') = a(a, a')q{a, a') 



1- / a{a,a')q{a,a')da' 



5(cj'-cj). (2) 



Convergence and ergodicity properties of the chain {X„} depend on the proposal 
kernel q, and they are studied extensively in \3T\. Jifc satisfies the DB condition with 
/ ensuring that it has stationary measure yi. is irreducible and aperiodic |32|, 
nonnegative definite, and reversible, thus the Markov chain with transition kernel 
converges in distribution to /i. 



2.1 Mixing times and speed of convergence 

It is known (7) that for a discrete-time Markov chain X„ with the transition kernel 
and the stationary distribution /, the rate of convergence to its stationarity can 
be measured in terms of kernel's second largest eigenvalue, according to 
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2||^"(a,.)-/l|rv<^^r 

where j3 ~ max{ |)3,„,„ | , j3i } and — 1 < )3„,„, <■■■ < )3o = 1 are the real eigenval- 
ues of ,Jiif . The spectral gap of kernel is defined by A — min | ; Var (h) ^ 
which for the self-adjoint kernel J(f , because of reversibility, is A(^) = 1 — j3i. 
With the Dirichlet form S and the variance defined by 

S{Kh) = \y. \Ko)-Ka'tJ€(^a,a')f{a), 
Var(/.) = i^|/.(a)-Ma')|V(cT')/(c^)- 

a, a' 

Between two algorithms producing Markov chains with identical equilibrium dis- 
tributions better in terms of the speed of convergence is the one with the smaller 
second eigenvalue in absolute value or equivalently with the larger spectral gap. 



3 The Coupled CGMC method 

The proposed algorithm is designed to generate samples from the microscopic prob- 
ability measure jj. with density / on a space E, coupling properly states of the mi- 
croscopic space E with states on a coarse space E having less degrees of freedom. A 
properly constructed coarse measure on E will be the basis for constructing efficient 
proposal kernels for MH algorithms sampling large systems. 

The coarsening procedure is based on the expansion of the target measure /i to a 
coarse and a finer part. Abstractly we write /(cr) = /(t] , ^ ) and E = E x E', where 
rj G E represents the coarse variables. 

We denote the projection operator on the coarse variables 

T:E^E, Ta = ri. 

The exact coarse marginal is 

Je' 

To obtain an explicit formula of the coarse marginal is as difficult as sampling the 
original target distribution since space £' remains high dimensional. Therefore use 
of approximating distributions of / becomes necessary. Such approximations have 
been proposed in | [T7][2Tj for stochastic lattice systems and are abstractly described 
in Section Q and for complex macromolecular systems see |4, 1 1 13||35). 

Denote /o an approximation of f on E. This distribution, combined with a re- 
construction distribution fr{£,\ri) corresponding to the finer variables ^, will con- 
struct a candidate for proposal distribution in MH algorithms performed in order 
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to sample from / at the original space E. An example of a 'good' proposal dis- 
tribution is /o((7) :— /o(tj)/, (^ I?])- For notational simplicity we write /r(<7|77) in- 
stead of /r(^|77). In terms of the Metropolis-Hastings algorithm this means that 
q{o,(j') = /o(c5''), or that /o(ct') is the stationary measure of the proposal kernel 
qia,a'). 

The coupled CGMC algorithm is composed of two coupled Metropolis iterations, 
the first generating samples from the proposal distribution and the second samples 
from the target measure. The first Metropolis step samples the coarse approximating 
marginal /o(i7), using an arbitrary proposal transition kernel ^o(j?,J?') to produce 
trial samples T]'. The second step is performed if the coarse trial sample is accepted, 
and consists of the reconstruction from the coarse trial state and a Metropolis crite- 
rion designed to ensure sampling from the correct microscopic density /. If a trial 
coarse sample is rejected, then we go back to the first step to rebuild a new coarse 
trial, so that the fine Metropolis step is not performed and no computational time is 
wasted on checking fine trial samples that are most likely to be rejected. 

In |j9) Efendiev et.al., propose the Preconditioning MCMC, a two stage ( coarse 
and fine ) Metropolis MCMC method, applied to inverse problems of subsurface 
characterization. The coarse and fine models are finite volume schemes of different 
resolutions for a PDE two-phase flow model. Our algorithm shares the same idea 
and structure with the Preconditioning MCMC of constructing a proposal density 
based on meso/macro-scopic properties of the model studied and taking advantage 
of the first stage rejections. In terms of the MC method 'coarsening' corresponds to 
enriching the range of the sampling measure based on coarse-scale models proposed 
by multiscale finite volume methods. The major difference of the Preconditioning 
MCMC and the proposed algorithm is that the latter alternates between different 
state spaces during each MC iteration, the coarse and the finer, whether in the for- 
mer the state space remains the same since coarse and fine problems are solved 
independently. Thus, at the end of a simulation we will have both fine-scale and 
"compressed", coarse-grained data. The performance of the coarse proposals in our 



case can be further estimated based on a systematic error analysis such as ( 14 1. 

The proposed procedure has also some common features with the modified Con- 
figurational bias Monte Carlo (CBMS) where the trial density is built up sequentialy 
with stage-wise rejection decision described in [27 1, applied effectively in quan- 
tum mechanical systems |5|. There are also some similarities with simulated sin- 
tering and transdimensional MCMC, see p7| and references therein. However, in 
our method, the construction of the variable dimensionality (and level of coarse- 
graining) state spaces and the corresponding Gibbs measures relies on statistical 
mechanics tools that allow a systematic control of the error from one level of coarse- 
graining to the next, e.g. ( [T4) l. 



3.1 The algorithm 



Coupled coarse graining and Markov Chain Monte Carlo for lattice systems 7 

We describe in detail the coupled CGMC Metropolis algorithm outlined in the 
previews section. 

Algorithm 2 ( Coupled CGMC Algorithm) 

Let Xq = (7o arbitrary, for n = 0, 1 , 2, . . . 
Given X„ = a 

Step 1 Compute the coarse variable r\ = Ta 
Step 2 Generate a coarse sample t/ ' ^ > T ' ) 
Step 3 Coarse Level Accept-Reject 
Accept Tj' with probability: 



If Tj' is accepted then proceed to Step 4 
else generate a new coarse sample Step 2 

Step 4 Reconstruct o' given the coarse trial r\'. 



Step 5 Fine Level Accept-Reject 
Accept o' with probability 



Steps 2 and 3 generate a Markov chain {Z„} in the coarse space £ with the 
transition kernel 



^('7,'?') = acG(T?,Tj')^o(T?,Tj') + 



acG(T?,z)^o(T?,z) 



5(T?'-TJ). 



The stationary measure of kernel ^ is fo{y\). Combination of this kernel and Steps 
1 and 4 constructs the desired proposal transition kemel qt^io^o'^ on the fine level 
space JL, 

qo{a,a') = ^{n,r]')fr{a'\r]'). 
According to the MH algorithm in order to sample from /, the fine level acceptance 
probabiUty should be af{a,a') = min |l, /(CT)^^p(a o^) }' ^^^^ ^ satisfies the 
Detailed Balance condition ^{r\,r\')f[^(r\) = =2(t/',t7)/o(t/'), a/ is equal to 



CCfio, o I = rmn < 1, — — ,. „ . — -. — - 



= min < 1 



) 

f{(y')Un)fr{o\n) 
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The chain X„ produced by Coupled CGMC algorithm is a Markov chain on the 
fine space Z, with the transition kernel 

5{a'-a). 

(3) 

The Markov chain X„ generated by the Coupled CGMC algorithm converges to 
the correct stationary distribution / and is ergodic, which ensures that ^ YJj=i hQ^j) 
is a convergent approximation of the averages f h{a)f{G)da for any h G L}{f)- 
Ergodicity and reversibility properties are satisfied ensuring that the algorithm gen- 
erates samples from the correct measure. We state this fact as a separate theorem 
proof of which is given in detail in |[15]. 

We denote E = {a € I.;f{o) > 0}, £ = {77 e t■JQ{^) > 0}. 

Theorem 1 For every conditional distribution qo, and fr such that the support of 
qofr includes E, 

i) The transition kernel satisfies the detailed balance (DB) condition with f 

jeccicy, o')f{a) = JTcgIct', o)f{a') 

ii) f is a stationary distribution of the chain. 

Hi) if qo{(J, <j') > 0, V(7, o' G E and E C supp^fo) then X„ is f -irreducible 
iv) is aperiodic 



1- / af{a,a')qoi(y,a')da' 



3.2 The rate of convergence 

The calculation of the rate of convergence to stationarity is a hard problem since it 
is model dependent as argued earlier What we can prove though for the proposed 
method is that it is comparable to the classical Metropolis-Hastings algorithm de- 
scribed in Algorithm[T] This fact is stated rigorously in the following theorem which 
we prove in [ 15 |. 

Let A(J^cg)j'^("^) be the spectral gap corresponding to the coupled CGMC 
J^CG, and the classical MH J^^-, (j2|, transition kernels respectively. 

Theorem 2 Let q{cy,cy') be a symmetric proposal transition probability for the 
classical MH algorithm and qo{Tl,ri') a symmetric proposal transition probability 
on the coarse space Efor the coupled CGMC algorithm, then for any reconstruction 
conditional probability fr{(y\r\) 

i) 

jeccicT, ct') = s^{(y, a')d§{a, a')Xc{a, a') (4) 
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^(CT,(7') 



q(a,&) 
M ■r\)fr(a\r\) 



Furthermore we define the subsets 



Ci 

C3 
C4: 



{{a,a') e E X E : {a < l.OcG < ^,<^f < ^} or {a = l,acG = l,"/ = l}} 
{{a,a') e E X E : {a = l,acG < '^^(^f = or {a < l,acG = < l}} 
{{a,a') e E X E : {a ^ l,acG =^ '^,ccf < or {a < l,acG < l,"/ = l}} 

{{a,a') e E X E : {a < l,acG ^ '^,0Cf ^ 1} or {a = l,acG < < l}} 



(a, a 



1, 



min{#^,#2i}, 
'-/o('7) ' Mn'l'' 



'/(cT,cj')eCi 
i/((7,(j')eC2 



/o('7) '/oM^, ,^ 
/(a)/o(T]') ' /(ct')/o(t?) J'' ^ ) ^3 



w/iere = inffj ^j/ (a, a') and 7 > 0, 7 > such that 7 < ^(c7, a') ^ 7- 



(5) 



4 Extended Lattice Systems 

This class of stochastic processes is employed in the modeling of adsorption, des- 
orption, reaction and diffusion of chemical species in numerous applied science 
areas such as catalysis, microporous materials, biological systems, etc. p][26l. To 
demonstrate the basic ideas, we consider an Ising-type system on a periodic d- 
dimensional lattice Aj^ with — n'' lattice points. At each x £ we can de- 
fine an order parameter a(x); for instance, when taking values and 1, it can de- 
scribe vacant and occupied sites. The energy Hn of the system, at the configuration 
O = {g{x) : X £ An} is given by the Hamiltonian, 

HN{cy) = -\ E J^[K{x-y)+J{x-y)]a{x)a{y)+Y,ha{x), (6) 

where h, is the external field and J is the inter-particle potential. Equilibrium states at 
the temperature ^ )3^' are described by the (canonical) Gibbs probability measure, 
and Za„ is the normalizing constant (partition function) 



Majv,/3 {d<y) = exp ( - pHNia))PN{da) . 



(7) 
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Furthermore, the product BernoulU distribution PAf(c7) is the prior distribution on 
An- 

The inter-particle potentials J account for interactions between occupied sites. 
We consider K corresponding to the short and J to the long-range interactions dis- 
cussed in detail in Section (4.2 1. General potentials with combined short and long- 
range interactions are discussed here, while we can also address potentials with 
suitable decay/growth conditions [2 J. 

The prior fV(iia) is typically a product measure, describing the system at 
/3 = 0, when interactions in are unimportant and thermal fluctuations-disorder- 
associated with the product structure of PN{da) dominates. By contrast at zero 
temperature, j5 = °o interactions, and hence order, prevail. Finite temperatures, 
< j3 < oo, describe intermediate states, including possible phase transitions be- 
tween ordered and disordered states. For both on-lattice or off-lattice particle sys- 
tems, the finite-volume equilibrium states of the system have the structure (|7]). 



4.1 Coarse-graining of microscopic systems 

Coarse-graining of microscopic systems is essentially an approximation theory and 
a numerical analysis question. However, the presence of stochastic fluctuations on 
one hand, and the extensivity of the models (the system size scales with the number 
of particles) on the other, create a new set of challenges. We discuss all these issues 
next, in a general setting that applies to both on-lattice and off-lattice systems. 

First, we write the microscopic configuration a in terms of coarse variables Tj and 
corresponding fine ones ^ so that a — [i] .^). We denote by T the coarse-graining 
map To = TJ. 

The CG system size is denoted by M, while the microscopic system size is 
= Mq, where we refer to q as the level of coarse graining, and q=\ corresponds 
to no coarse graining. The exact CG Gibbs measure is given (with a slight abuse of 
notation) by /iM,/3 = IJ-n,^ o r^' . In order to write IXm.P ^ more convenient form 
we first define the CG prior pMidrj) ^ Pn° T^^. The conditional prior Pf^{dG\ri) is 
the probability of having a microscopic configuration a, given a coarse configura- 
tion TJ . We now rewrite p.^ p using the exact coarse-grained Hamiltonian: 

^-PHmM ^ E[e-I^"''\r]]= J e-P"''^"^PNida\ri), (8) 

a procedure known as the renormalization group map, |jT2j; flM.pidri) is now re- 
written using (|8]l as 



flM,p{dll) = ^e-^"-^^^PM{dll). 



(9) 
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Although typically PM{dr\) is easy to calculate, even for moderately small values 
of the exact computation of the coarse-grained Hamiltonian 5m(^7) given by (|8]l 
is, in general, impossible. 



We have shown in |21 1 that there is an expansion of HM{ri) into a convergent 
series 

Hm{i) = Hj°\ii)+Hj^\l) +Hj^\n) + ■■■ +Hl^\n)+N X ^(e'^) , (10) 

by constructing a suitable first approximation Hj^\ri) and identifying a suitable 
small parameter e to control the higher order terms in the expansions. Truncations 
including the first terms in ([TO]l correspond to coarse-graining schemes of increasing 
accuracy. In order to obtain this expansion we rewrite ([8| as 



//m(t7) -//^(r?) -^\ogE[e-P^"--"M^^^^\ri] . (11) 



We need to show that the logarithm can be expanded into a convergent series, yield- 
ing eventually ( [TOj i, however, two interrelated difficulties emerge immediately: first, 
the stochasticity of the system in the finite temperature case, yields the nonlin- 
ear log expression which in turn will need to be expanded into a series. Second, 
the extensivity of the microscopic system, i.e., typically the Hamiltonian scales as 
Hn = 0{N), does not allow the expansion of the logarithm and exponential func- 
tions into a Taylor series. For these two reasons, one of the mathematical tools 
we employed is the cluster expansion method, see p3| for an overview. Clus- 
ter expansions allow us to identify uncorrelated components in the expected value 

E[g-/3(^fA'-WM |tj] , which in turn will permit us to factorize it, and subsequently 
expand the logarithm. 



The coarse-graining of systems with purely long- or intermediate-range interac- 
tions of the form 

J{x-y)=L-^v[{x-y)/Ly x,yeAN, (12) 

where V{r) = V{—r), V(r) ~ 0,\r\ > 1, was studied using cluster expansions in 
|j2^^20 21 1. The corresponding CG Hamiltonian is 



(13) 



leAM ''<^^M leAM keA„ 



j{k,i) = ^Y.L J(^-y)^ - L ■^(^-3') 
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One of the results therein is on deriving error estimates in terms of the specific 
relative entropy ^{lJ.\v) := N^^Y.a^og{ii{a)/v{a)}lJ.{a) between the corre- 
sponding equilibrium Gibbs measures. Note that the scaling factor A^^' is related to 
the extensivity of the system, hence the proper error quantity that needs to be tracked 
is the loss of information per particle. Using this idea we can assess the information 
compression for the same level of coarse graining in schemes differentiated by the 
truncation level p in ( [T0| 

mi^^p\^N.poT-')^ff{e''+'), e^^||Vy||i(|), (14) 

where Hj^\ri) in ([TO]i is given by ( [T3] ). The role of such higher order schemes was 
demonstrated in nucleation, metastability and the resulting switching times between 
phases, |2|. 

Although CGMC and other CG methods can provide a powerful computational 
tool in molecular simulations, it has been observed that in some regimes, important 
macroscopic properties may not be captured properly. For instance, (over-)coarse 
graining in polymer systems may yield wrong predictions in the melt structure | ll; 
similarly wrong predictions on crystallization were also observed in the CG of com- 
plex fluids, |31 1. In CGMC for lattice systems, hysteresis and critical behavior may 
also not be captured properly for short and intermediate range potentials, p9|[2T[ . 
Motivated by such observations, in our recent work we studied when CG methods 
perform satisfactorily, and how to quantify the CG approximations from a numerical 
analysis perspective, where error is assessed in view of a specified tolerance. Next, 



we discuss systems with long range interactions, i.e., L >> 1 in ( 12 1. These systems 
can exhibit complex behavior such as phase transitions, nucleation, etc., however, 
they are more tractable analytically. At the same time they pose a serious challenge 
to conventional MC methods due to the large number of neighbors involved in each 
MC step. 

Here we adopt this general approach, however, the challenges when both short 
and long-range interactions are present, require a new methodology. Short-range 
interactions induce strong "sub-coarse grid" fine-scale correlations between coarse 
cells, and need to be explicitly included in the initial approximation h]^'(tj). For 



this reason we introduced in p4) a multi-scale decomposition of the Gibbs state 
(|7]), into fine and coarse variables, which in turn allows us to describe in an explicit 
manner the communication across scales, for both short and long-range interactions. 



4.2 Multiscale Decomposition and Splitting Methods for MCMC 

We first focus on general lattice systems, and subsequently discuss related applica- 
tions in later sections. We consider (j6|l where in addition to the long-range potential 
([T2|, we add the short-range K{x~y) = S'^U {N\x-y\/S), where S « L and U 
has similar properties as V in ( [T2] i; for 5 = 1 we have the usual nearest neighbor 
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interaction. The new Hamiltonian includes both long and short-range interactions: 

The common theme is the observation that long-range interactions L >> 1 can 
be handled very efficiently by CGMC, ([T4]i. On the other hand short-range interac- 
tions are relatively inexpensive and one could simulate them with Direct Numerical 
Simulation (DNS) provided there is a suitable splitting of the algorithm in short 
and long-range parts, that can reproduce within a given tolerance equilibrium Gibbs 
states and dynamics. We return to the general discussion in ( fTO] ) and outline the 
steps we need in order to construct the CG Hamiltonian for the combined short and 
long-range interactions. 

Step 1: Semi-analytical splitting schemes. Here we take advantage of CG approxi- 
mations developed in ([T4]i in order to decompose our calculation into analytical and 
numerical components, the latter involving only short-range interactions: 



^iN,p{da)^e-l^"^^''^PNida)^ 



where h'j^ is the analytical CG formula ( [T3] l constructed for the computationally 
expensive, for conventional MC, long-range part; due to the estimates ( [l4j l, the 
first term has controlled error Furthermore, the dependence of e on W in these 
estimates suggests a rearrangement of the overall combined short- and long-range 
potential, into a new short-range interaction that includes possible singularities orig- 
inally in the long-range component ( [T2] i, e.g., the singular part in a Lennard- Jones 
potential, and a locally integrable (or smooth) long-range decaying component that 
can be analytically coarse-grained using ( [T3] l, with a small error due to ( [l4j i. This 
breakdown allows us to isolate the short-range interactions (after a possible re- 
arrangement!), and suggests the two alternative computational approaches: either 
seek an approximation c^^^mM — J e^P^NPi^(^d(y\ri), or use sampling methods to 
account for the short-range "unresolved" terms. 



4.3 Microscopic Reconstruction 

The reverse procedure of coarse-graining, i.e., reproducing "atomistic" properties, 
directly from CG simulations is an issue that arises extensively in the polymer sci- 
ence literature, 1 30 36). The principal idea is that computationally inexpensive CG 



simulations will reproduce the large scale structure and subsequently microscopic 
information will be added through microscopic reconstruction, e.g., the calculation 
of diffusion of penetrants through polymer melts, reconstructed from CG simula- 
tion, f30l. In this direction, CGMC provides a simpler lattice framework to math- 
ematically formulate microscopic reconstruction and study related numerical and 
computational issues. Interestingly this issue arised also in the mathematical error 
analysis in 1 18 22) . 
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The mathematical formulation for the reconstruction of the microscopic equilib- 
rium follows trivially when we rewrite the Gibbs measure (|7]) in terms of the exact 
CG measure corresponding to (|8]l, defined in (|9|, |20j : 

We can define the conditional probabiUty iJ.N{d(j\ri) as the exact reconstruction 
of jiNida) from the exactly CG measure fLMidrj). Although many fine-scale con- 
figurations a correspond to a single CG configuration Tj, the "reconstructed" con- 
ditional probability iiN{da\ri) is uniquely defined, given the microscopic and the 
coarse-grained measures ^^{da) and fLMidrj) respectively. 

A coarse-graining scheme provides an approximation p.^^{dri) for fiMidrj), at 
the coarse level. The approximation il^^{dri) could be, for instance, any of the 
schemes discussed in Section |4.2| To provide a reconstruction we need to lift the 
measure jl^'^{dri) to a measure fi^'^{da) on the microscopic configurations. That 
is, we need to specify a conditional probability VN{dG\ri) and set n^'^{dG) :~ 
VN{dG\ri)jl^^{dri) . In the spirit of our earlier discussion, it is natural to measure 
the efficiency of the reconstruction by the relative entropy, 

^(^riM ='^(A;''lAM)+/^(v^(-|J7)IM-|T7))Ar(«'^7), (15) 

i.e., relative entropy spUts the total error at the microscopic level into the sum of the 
error at the coarse level and the error made during reconstruction, |20 34 1. 

The first term in ( fTSj ) can be controlled via CG estimates, e.g., ( [I4j l. However, 
(T5\ suggests that in order to obtain a successful reconstruction we then need to 
construct VN{da | ?]) such that (a)^{vNida \ rj) \^iN{da \ 7])) should be of the same 
order as the first term in ([T5|, and (b) it is easily computable and implementable. 

The simplest example of reconstruction is obtained by considering a microscopic 



system with intermediate/long-range interactions ( 12 1 



fi'J^{dri)^fLj^\dii), VN{da\n)^PN{da\ii). (16) 

Thus we first sample the CG variables 7] involved in using a CGMC algorithm; 
then we reconstruct the microscopic configuration a by distributing the particles 
uniformly on the coarse cell, conditioned on the value of Tj. Since P!^{da\ri) is a 
product measure this can be done numerically in a very easy way, without com- 
munication between coarse cells and only at the coarse cells where an update has 
occurred in the CGMC algorithm. In this case the analysis in [ ,23) yields the esti- 
mates 

^(Mri/iM)=0(e2), ^(^^(.|T,)|P^(.|Tj)) = ^(//(«)(77)-//(77)) =0(e2). 



Hence the reconstruction is second order accurate and of the same order as the 
coarse-graining given by ( [T3| l. 
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5 Example: Short and long-range interactions 



Short and long-range interactions pose a formidable computational challenge. We 
consider an example that has been explicitly solved by Kardar in |16|. The model 
considered has state space Z^f — 

{0,1}'^'^^ where Aa, is a 1 -dimensional lattice with 
sites. The energy of the system at configuration a = {a{x),x G A^f} is 



^ X \x-y\ = \ X y^x 

= H-;,{a)+H'^{a)+E{a). 

Hamiltonian H^io) consists of the short-range term Hj^, the long-range term //^ 
and an external field E. The interactions involved in are of the nearest-neighbor 
type with strength while //j^ represents a mean-field approximation or the Curie- 
Weiss model defined by the potential J averaged over all lattice sites. For this generic 
model Kardar gave in [ 16j a closed form solution for magnetization {K,J,h),foi 
the state space {—1,1} 



Mj3 {K,J,h) — cLTg min [^rr? — log 



e^cosh(/! +/m) + W e^^ sin^ {h + Jm) + e^^^ 



a simple rescaling of which gives the exact average coverage mp{K,J,h) for the 
lattice-gas model considered here. 

mp{K,J,h) = i (mp (^^K, iy, \h-^-J- ^K^ + 1^ (17) 

We have constructed the classical single spin-flip M-H algorithm and the coupled 
Metropolis CGMC for the single spin-flip algorithm, both generating samples from 
the Gibbs measure 



Z, 



'N 



We denote G'^ the state that differs from a only at the site x, G-'^{y) = G{y),y ^ 
X, (7^{x) 1 — cj(x), the proposal transition kernel is q{a'\<j) = ' ~ f^)' 

proposing a spin-flip at the site x with the probability ^ . 

We apply the coupled CGMC method with coarse updating variable 

ri:=T(j = {ri{k),k^l,...,M} 



^(^) •= I^rGQ = N with a coarsening level q < M, and for the maximum 

coarsening q=N where the coarse variable is total magnetization r\ — Y^xeA <^{x)- 
This can be thought as a coarsening procedure constructing a system consisting of 
one big coarse cell M — I with q = N sites. Since we want to consider only single 
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spin-flip updates, for the sake of comparison to the classical Metropolis method, the 
cell updating can take only the values ±1 and the reconstruction is chosen uniform 



in each cell, in the sence described in example at Section 4.3 



Table[T]gives a comparison of the classical single-site updating Metropolis Hast- 
ings algorithm with the proposed coupled Metropolis CGMC algorithm, in terms 
of computational complexity per iteration. By computational complexity here we 
mean the cost of calculating energy differences involved at the acceptance proba- 
bilities. Consider the case that both the microscopic single-site updating Metropo- 
lis and the two-step CGMC are run n times. This is reasonable to consider since 
as stated at Theorem |2] the two methods have comparable mixing times, therefore 
the number of iterations needed to achieve stationarity are comparable. We denote 
'^{(k:G) '■— 1 1 0!cg(J?, J?')^o('?' ^')/o('?)'^'?'^^' the average acceptance rate of the 
coarse proposal. The average number of accepted coarse samples is «i :~ [E(acG)«] 
for which «i < n since E(acG) < 1 • This means that the reconstruction and fine step 
acceptance criterion are performed in average only for rii iterations. 



Table 1 Operations count for evaluating energy differences for n iterations 

Cost Metropolis Hastings Coupled CGMC q <N Coupled CGMC q = N 

coarse A-R - nxO{M) nxO{l) 

fineA-R nxO{N) ni x 0(1) ni x 0(1) 




Fig. 1 Phase Diagram 







Results of computational implementation are shown in Figure |2] and Table |2] and 
[3] Figure represents the average coverage versus the external field h for the exact 
solution mj.v, the classical MH result < nid > and the coupled CGMC < m >, for 
a choice of interaction parameters K = 1, 7 = 5 in the ferromagnetic region as 
is stated at the phase diagram depicted in Figure [T| The exact solution m^^ as is 
ploted in Figure corresponds to the part the full solution (17 1 up to the point it 
jumps. Figure |2V) is a graph of the average acceptance rates for the classical MH 
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Fig. 2 N= 1028, q=8, K= 1, J= 5: (a) Coverage ; (b) Average acceptance. 




-6 -4 -2 2 4 6 8 10 12 
External field h 



Fig. 3 N= 1028, q=8, K= 1, J= 5: Local error log(| <m> - < m,^\) 

algorithm and the coupled CGMC algorithm, that verifies the theoretical proof of 
the fact that the two algorithms have comparable mixing times since the acceptance 
rate is strongly related to mixing times. In the same figure we also give the average 
acceptance rates of the coarse and fine step of the coupled method, noting that the 
fine acceptance rate is high which means that most of the trial samples entering the 
fine step are accepted. 

Table [2] reports the error between the exact solution and the average cover- 
age obtained from the coupled CGMC algorithm. Error is measured in terms of 

the pointwise solutions as Error^. — {Y.i{'^yx{hi)— < m > {hi) f) and Error,., = 

ij^i{mex[hi)— < m^-i > (/z,))^) '''^ for the coupled and the classical method respec- 
tively, where hi are the different external field parameters for which the average 
coverages are computed. CPU times are compared for the coarse-graining levels 
q = A and ^ = 8. To demonstrate the robustness of the algorithm we present simu- 
lations at three different points of the phase diagram plane K J: in the disordered 
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Table 2 N=4096 





CG level q 


Errorc 


CPU(min) 


K= -2.0, 3=2 


4 
8 


0.089 
0.302 


93.52 
45.8 


K= 1.0, J=5 


4 
8 


0.003 
0.003 


93.6 
45.9 


K= 1, J=l 


4 


0.027 


91.6 


8 


0.100 


45.5 



((K = -2.0,7 ^ 2) and = 1,7 = 1)) and ferromagnetic (K = 1.0,7 = 5) regions. 
In table [3] we compare the error between the coupled CGMC average coverage with 
the exact solution and the corresponding CPU time for ^ = 4 and ^ = 8, in the 
ferromagnetic region (K = 1.0,7 = 5) for which the classical Metropolis results. 

These results demonstrate the efficiency of the coupled CGMC methods in terms 
of computational time since the run time gain scales almost linearly with the coars- 
ening level. We expect that according to Theorem[2|i, the error should be indepen- 
dent of the coarse graining parameter due to the microscopic nature of the algo- 
rithm though this is not evident in the tables since we are using a simplification of 
the reconstruction procedure for computational ease. We should also mention that a 
large number of samples (10^) were considered ensuring the statistical error is small 
enough. 



Table 3 N= 1028, K= 1, 1= 5, Error,.; = 0.003, Cla.ssical CPU = 94.5min 



CG level Error,. Coupled CPU(min) 



q=4 0.01 23.1 
q=8 0.04 12.1 



6 Conclusions 

An advantage of the Coupled CGMC approach over the asymptotics methodology 



discussed in Section 4.2 is that the trial distribution may even be order one away 
from the target distribution, however, the method can still perform well. On the other 
hand, the methods can complement each other; for example, for equilibrium sam- 
pUng considered in this work we use as a trial reconstructed distribution, the con- 
ditional measure v{da\r\) in the multiscale decomposition in |24|, see also Section 



4.3 Such proposals based on careful statistical mechanics-based approximations 
provide better trial choices for the MH methods and more efficient sampling, as is 
proved theroretically and numerically. The example illustrated makes clear that the 
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coupled CGMC method implements a splitting of the short and long-range interac- 
tion terms, into the two Metropolis acceptance criteria involved. The long-range part 
which is responsible for the expensive calculations at a fully microscopic method, 
now enters only in the coarse approximation measure where its computational cost 
is much lower. 

Coupling of a coarse and fine step is also effective in the study of dynamic pro- 
cesses of stochastic lattice systems with kinetic Monte Carlo methods, a topic stud- 
ied in detail in fT5\ . 
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